% IHMS GHR 2024.01.08	该程序实现 FBCCA 算法，在主程序中需要初始化 main.m 中的参数。

function [outputCommand, correlationCoefficient, correlationCoefficient1, correlationCoefficient2, costNow] = FBCCA_algorithm(EEGDataNow, ...
    ReferenceSignal, ChannelNumber, DownSampleFactor, CorrectionFactor, FrequencySample, WindowDownsample, StimulationFrequency, ...
    CommandList, FilterBankFactor, CCFlag, WeightedMeanFactor, Threshold)
    analysisData = EEGDataNow(1:ChannelNumber,:);
    % 《《==================== 预处理 ====================》》
    dataDownsample = downsample(analysisData', DownSampleFactor)';  % 降采样
    dataCorrection = BaselineCorrection(CorrectionFactor, ChannelNumber, WindowDownsample, dataDownsample);  % 基线修正
    dataBandpass = BandpassFilter(FrequencySample, dataCorrection);  % 带通滤波
    dataBandstop = BandstopFilter(FrequencySample, dataBandpass);  % 带阻滤波

    % 《《==================== 数据处理 ====================》》
    dataFilterBank = FilterBank(FilterBankFactor, FrequencySample, dataBandstop, StimulationFrequency);  % 滤波器组滤波
    [outputCommand, correlationCoefficient, correlationCoefficient1, correlationCoefficient2, costNow] = CCAAnalysis(FilterBankFactor, ...
        WeightedMeanFactor, dataFilterBank, ReferenceSignal, CommandList, CCFlag, Threshold);  % 典型相关分析
end


% 《《==================== 函数体 ====================》》
% 基线修正 左右修正 *
function datacorrection = BaselineCorrection(n, nchan, windown, data)
    datacorrection = zeros(nchan, windown);
    for i = 1:nchan
    average1 = mean(data(i, 1:n));
    average2 = mean(data(i, (windown-n):windown));
        for j = 1:windown
            if j <= n/2
                datacorrection(i, j) = data(i, j) - average1;
            elseif j > n/2 && j <= (windown-n/2)
                datacorrection(i, j) = data(i, j) - mean(data(i, (j-n/2):(j+n/2)));
            else
                datacorrection(i, j) = data(i, j) - average2;
            end
        end
    end
end

% 带通滤波 8~88 Hz *
function databandpass = BandpassFilter(fs, data)
	fp1 = [8, 88];  fs1 = [4, 100];
    Fs2 = fs/2;
    Wp = fp1/Fs2;  Ws = fs1/Fs2;
    Rp = 1;  Rs = 30;
    [n, Wn] = cheb2ord(Wp, Ws, Rp, Rs);
    [b1, a1] = cheby2(n, Rs, Wn);
    databandpass = filter(b1, a1, data, [], 2);
end

% 带阻滤波 50 Hz *
function databandstop = BandstopFilter(fs, data)
    fp1 = [4, 100];  fs1 = [49, 51];
    Fs2 = fs/2;
    Wp = fp1/Fs2;  Ws = fs1/Fs2;
    Rp = 1;  Rs = 2;
    [n, Wn] = cheb2ord(Wp, Ws, Rp, Rs);
    [b2, a2] = cheby2(n, Rs, Wn, 'stop');
    databandstop = filter(b2, a2, data, [], 2);
end

% 滤波器组带通滤波 *
function datafilterbank = FilterBank(N, fs, data, fk)
    datafilterbank = zeros(N, size(data, 1), size(data, 2));  % 分配储存空间
    for i = 1:N
        fp1 = [i*min(fk), 88];  fs1 = [i*min(fk)-i*1, 100];
        Fs2 = fs/2;
        Rp = 1;  Rs = 35;  % Rp 和 Rs 分别代表通带起伏和阻带衰减
        Wp = fp1/Fs2;  Ws = fs1/Fs2;
        [n, Wn] = cheb2ord(Wp, Ws, Rp, Rs);
        [b1, a1] = cheby2(n, Rs, Wn);
        datafilterbank(i, :, :) = filter(b1, a1, data, [], 2);
    end
end

% 典型相关分析处理数据（用 MATLAB 官方提供的函数 ） *
function [outputcommand, correlationcoefficient, correlationcoefficient1, correlationcoefficient2, costNow] ...
    = CCAAnalysis(N, ab, data, ref, cmd, flag, th)
    correlationcoefficient1 = zeros(size(cmd, 2), N);
    correlationcoefficient2 = zeros(size(cmd, 2), 1);
    for i = 1:size(cmd, 2)
        for j = 1:N
            [~,~,r,~,~] = canoncorr(squeeze(data(j, :, :))', squeeze(ref(i, :, :))');
            correlationcoefficient1(i, j) = r(1);
        end
        if flag == 0
            correlationcoefficient2(i) = norm(squeeze(correlationcoefficient1(i, :)));  % 取模
        else
            correlationcoefficient2(i) = WeightedMean(N, ab(1), ab(2), squeeze(correlationcoefficient1(i, :)));  % 加权平均
        end
    end
	[correlationcoefficient, I] = max(correlationcoefficient2);
	outputcommand = cmd(I);
    costNow = ResultCredibility(correlationcoefficient2, I, th);  % 计算当前命令的代价值
end

% 加权平均 *
function result = WeightedMean(n, a, b, data)
    result = 0;
    for i = 1:n
        wn = (i^(-a) + b);
        corrcoe = wn * data(i)^2;
        result = result + corrcoe;
    end
end

% 计算当前命令的可靠程度 *
function output = ResultCredibility(cc, I, th)
    output = zeros(3, 1);
    % 获取第二大特征值的位置
    I1 = I;  cc1 = cc;
    cc1(I1) = nan;
    [~, I2] = max(cc1);
	Softcc = softmax(cc);  % 将相关系数使用 Softmax 函数处理
    % 计算代价值
    CostH0 = 0;  CostHq = 0;
    for i = 1:size(Softcc, 1)
        if i == I1
            C = 1;
        elseif i == I2
            C = -1;
        else
            C = 0;
        end
        CostH0 = CostH0 + th * log(Softcc(i));
        CostHq = CostHq + C * log(Softcc(i));
    end
    output(1) = -CostH0;  output(2) = -CostHq;
    if output(1) < output(2)
        output(3) = 1;
    end
end

